suppressMessages(require(Seurat))
suppressMessages(require(ggplot2))
suppressMessages(require(cowplot))
suppressMessages(require(scater))
suppressMessages(require(scran))
suppressMessages(require(BiocParallel))
suppressMessages(require(BiocNeighbors))
setwd("/home/li/Transcriptomic_patients/gse145926")
#draw the plot group.by celltype
hms_individual_integrated<-readRDS(file="hms_individual_integrated_OK.rds")
p1 <- DimPlot(hms_individual_integrated, reduction = "umap", group.by = "celltype")
p1
#find how many 15cluster
ElbowPlot(hms_individual_integrated)
hms_neighbor<- FindNeighbors(hms_individual_integrated, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
hms_cluster <- FindClusters( hms_neighbor, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 52158
## Number of edges: 1604563
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8838
## Number of communities: 15
## Elapsed time: 18 seconds
head(Idents(hms_cluster), 5)
## h1_AAACCTGAGACACTAA-1 h1_AAACCTGAGGAGTACC-1 h1_AAACCTGAGGATATAC-1
## 1 9 1
## h1_AAACCTGAGGTCATCT-1 h1_AAACCTGCACGGATAG-1
## 1 3
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
hms_cluster<- RunUMAP(hms_cluster, dims = 1:10)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 03:55:40 UMAP embedding parameters a = 0.9922 b = 1.112
## 03:55:40 Read 52158 rows and found 10 numeric columns
## 03:55:40 Using Annoy for neighbor search, n_neighbors = 30
## 03:55:40 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 03:55:46 Writing NN index file to temp file /tmp/Rtmp3WFr3g/file36f860cdd664
## 03:55:47 Searching Annoy index using 1 thread, search_k = 3000
## 03:56:05 Annoy recall = 100%
## 03:56:05 Commencing smooth kNN distance calibration using 1 thread
## 03:56:07 Initializing from normalized Laplacian + noise
## 03:56:09 Commencing optimization for 200 epochs, with 2146824 positive edges
## 03:56:28 Optimization finished
DimPlot(hms_cluster, reduction = "umap")
saveRDS(hms_cluster, file = "hms_cluster_test.rds")
#name each cluster id
new.cluster.ids <- c("Macrophage", "Macrophage", "Macrophage", "Macrophage", "Neutrophils", "Macrophage", "Naive_CD4_T", "NK", "Neutrophils", "Dendritic", "T", "T","Basal","Plasma","T")
names(new.cluster.ids) <- levels(hms_cluster)
hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)
DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.
saveRDS(hms_cluster_id, file = "hms_cluster_id_test.rds")
#hms_cluster_id<-readRDS(file="hms_cluster_id_test.rds")
Macrophage<-subset(hms_cluster_id, idents=c('Macrophage'))
DimPlot(Macrophage, reduction = "umap")
saveRDS(Macrophage, file="Macrophage.rds")
Neutrophils<-subset(hms_cluster_id, idents=c('Neutrophils'))
DimPlot(Neutrophils, reduction = "umap")
saveRDS(Neutrophils, file="Neutrophils.rds")
Naive_CD4_T<-subset(hms_cluster_id, idents=c('Naive_CD4_T'))
DimPlot(Naive_CD4_T, reduction = "umap")
saveRDS(Naive_CD4_T, file="Naive_CD4_T.rds")
NK<-subset(hms_cluster_id, idents=c('NK'))
DimPlot(NK, reduction = "umap")
saveRDS(NK, file="NK.rds")
Dendritic<-subset(hms_cluster_id, idents=c('Dendritic'))
DimPlot(Dendritic, reduction = "umap")
saveRDS(Dendritic, file="Dendritic.rds")
Basal<-subset(hms_cluster_id, idents=c('Basal'))
DimPlot(Basal, reduction = "umap")
saveRDS(Basal, file="Basal.rds")
T<-subset(hms_cluster_id, idents=c('T'))
DimPlot(T, reduction = "umap")
saveRDS(T, file="T.rds")
Plasma<-subset(hms_cluster_id, idents=c('Plasma'))
DimPlot(Plasma, reduction = "umap")
saveRDS(Plasma, file="Plasma.rds")
#input each cluster
Basal<-readRDS("Basal.rds")
Dendritic<-readRDS("Dendritic.rds")
Macrophage<-readRDS("Macrophage.rds")
Naive_CD4_T<-readRDS("Naive_CD4_T.rds")
Neutrophils<-readRDS("Neutrophils.rds")
NK<-readRDS("NK.rds")
Plasma<-readRDS("Plasma.rds")
T<-readRDS("T.rds")
DimPlot(Basal, reduction = "umap", split.by = "tech")
DimPlot(Dendritic, reduction = "umap", split.by = "tech")
DimPlot(Macrophage, reduction = "umap", split.by = "tech")
DimPlot(Naive_CD4_T, reduction = "umap", split.by = "tech")
DimPlot(Neutrophils, reduction = "umap", split.by = "tech")
DimPlot(NK, reduction = "umap", split.by = "tech")
DimPlot(Plasma, reduction = "umap", split.by = "tech")
DimPlot(T, reduction = "umap", split.by = "tech")
markers.to.plot <- c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","ARRDC3","EREG","ARSE","MSP2K6","DHCR7","UCP2","SLC25A10","VIL1","MCM5","DHCR24","SLC9A3R1","PFN1","TPPP3","DEGS2","RAB26")
DoHeatmap(subset(Basal,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(Basal, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(Dendritic,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(Dendritic, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(Macrophage,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(Macrophage, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(Naive_CD4_T,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(Naive_CD4_T, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(Neutrophils,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(Neutrophils, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(NK,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(NK, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(Plasma,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(Plasma, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
DoHeatmap(subset(T,downsample=50000), features = markers.to.plot, size = 5, group.by="tech")
## Warning in DoHeatmap(subset(T, downsample = 50000), features =
## markers.to.plot, : The following features were omitted as they were not found
## in the scale.data slot for the integrated assay: RAB26, DEGS2, PFN1, SLC9A3R1,
## MCM5, VIL1, SLC25A10, DHCR7, MSP2K6, ARSE, ARRDC3
RidgePlot(Basal, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"),cols = c("green3","cornflowerblue","orangered"), group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.00325
## Picking joint bandwidth of 0.0474
## Picking joint bandwidth of 0.117
## Picking joint bandwidth of 0.0865
## Picking joint bandwidth of 0.184
## Picking joint bandwidth of 0.17
## Picking joint bandwidth of 0.00184
## Picking joint bandwidth of 0.0186
## Picking joint bandwidth of 0.0215
## Picking joint bandwidth of 0.268
## Picking joint bandwidth of 0.0309
## Picking joint bandwidth of 0.153
## Picking joint bandwidth of 0.0523
## Picking joint bandwidth of 0.216
RidgePlot(Dendritic, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"),cols = c("green3","cornflowerblue","orangered"), group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.00878
## Picking joint bandwidth of 0.0276
## Picking joint bandwidth of 0.0504
## Picking joint bandwidth of 0.0327
## Picking joint bandwidth of 0.0744
## Picking joint bandwidth of 0.185
## Picking joint bandwidth of 0.00438
## Picking joint bandwidth of 0.0285
## Picking joint bandwidth of 0.0113
## Picking joint bandwidth of 0.22
## Picking joint bandwidth of 0.0226
## Picking joint bandwidth of 0.139
## Picking joint bandwidth of 0.00564
## Picking joint bandwidth of 0.0065
RidgePlot(Macrophage, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"), cols = c("green3","cornflowerblue","orangered"),group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.0032
## Picking joint bandwidth of 0.0547
## Picking joint bandwidth of 0.0595
## Picking joint bandwidth of 0.0232
## Picking joint bandwidth of 0.0672
## Picking joint bandwidth of 0.0975
## Picking joint bandwidth of 0.0025
## Picking joint bandwidth of 0.0244
## Picking joint bandwidth of 0.00684
## Picking joint bandwidth of 0.118
## Picking joint bandwidth of 0.027
## Picking joint bandwidth of 0.0945
## Picking joint bandwidth of 0.00494
## Picking joint bandwidth of 0.00208
RidgePlot(Naive_CD4_T, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"), cols = c("green3","cornflowerblue","orangered"),group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.000519
## Picking joint bandwidth of 0.007
## Picking joint bandwidth of 0.137
## Picking joint bandwidth of 0.0131
## Picking joint bandwidth of 0.105
## Picking joint bandwidth of 0.151
## Picking joint bandwidth of 0.00046
## Picking joint bandwidth of 0.00276
## Picking joint bandwidth of 0.0178
## Picking joint bandwidth of 0.156
## Picking joint bandwidth of 0.003
## Picking joint bandwidth of 0.094
## Picking joint bandwidth of 0.00795
## Picking joint bandwidth of 0.00122
RidgePlot(Neutrophils, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"), cols = c("green3","cornflowerblue","orangered"),group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.00741
## Picking joint bandwidth of 0.125
## Picking joint bandwidth of 0.174
## Picking joint bandwidth of 0.099
## Picking joint bandwidth of 0.208
## Picking joint bandwidth of 0.296
## Picking joint bandwidth of 0.00974
## Picking joint bandwidth of 0.145
## Picking joint bandwidth of 0.0618
## Picking joint bandwidth of 0.241
## Picking joint bandwidth of 0.0378
## Picking joint bandwidth of 0.0676
## Picking joint bandwidth of 0.00527
## Picking joint bandwidth of 0.00183
RidgePlot(NK, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"), cols = c("green3","cornflowerblue","orangered"),group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.000989
## Picking joint bandwidth of 0.00466
## Picking joint bandwidth of 0.127
## Picking joint bandwidth of 0.0269
## Picking joint bandwidth of 0.0874
## Picking joint bandwidth of 0.123
## Picking joint bandwidth of 0.00018
## Picking joint bandwidth of 0.00249
## Picking joint bandwidth of 0.00739
## Picking joint bandwidth of 0.144
## Picking joint bandwidth of 0.00219
## Picking joint bandwidth of 0.115
## Picking joint bandwidth of 0.0115
## Picking joint bandwidth of 0.00101
RidgePlot(Plasma, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"), cols = c("green3","cornflowerblue","orangered"),group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.00359
## Picking joint bandwidth of 0.016
## Picking joint bandwidth of 0.0848
## Picking joint bandwidth of 0.0223
## Picking joint bandwidth of 0.124
## Picking joint bandwidth of 0.171
## Picking joint bandwidth of 0.00139
## Picking joint bandwidth of 0.00849
## Picking joint bandwidth of 0.00974
## Picking joint bandwidth of 0.217
## Picking joint bandwidth of 0.00637
## Picking joint bandwidth of 0.162
## Picking joint bandwidth of 0.019
## Picking joint bandwidth of 0.0038
RidgePlot(T, features = c("IL1A", "CXCL2","TNFAIP3","MAFF","PPP1R15A","NFKBIA","PTX3","CXCL3","CCL20","IFIT2","EREG","UCP2","DHCR24","TPPP3"), cols = c("green3","cornflowerblue","orangered"),group.by="tech", ncol = 3) + theme(axis.title.y = element_blank())
## Picking joint bandwidth of 0.00406
## Picking joint bandwidth of 0.0337
## Picking joint bandwidth of 0.115
## Picking joint bandwidth of 0.036
## Picking joint bandwidth of 0.127
## Picking joint bandwidth of 0.159
## Picking joint bandwidth of 0.00745
## Picking joint bandwidth of 0.0229
## Picking joint bandwidth of 0.0128
## Picking joint bandwidth of 0.196
## Picking joint bandwidth of 0.0213
## Picking joint bandwidth of 0.139
## Picking joint bandwidth of 0.0111
## Picking joint bandwidth of 0.00282
